Accretion Disks Around Binary Black Holes: A Simple GR-Hybrid Evolution Model 
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We consider a geometrically thin, Keplerian disk in the orbital plane of a binary black hole 
(BHBH) consisting of a spinning primary and low-mass secondary (mass ratio g ^ 1). To account 
for the principle effects of general relativity (GR), we propose a modification of the standard Newto- 
nian evolution equation for the (orbit-averaged) time-varying disk surface density. In our modified 
equation the viscous torque in the disk is treated in full GR, while the tidal torque is handled in the 
Newtonian limit. Our GR-hybrid treatment is reasonable because the tidal torque is concentrated 
near the orbital radius of the secondary and is most important prior to binary-disk decoupling, when 
the orbital separation is large and resides in the weak-field regime. The tidal torque on the disk 
diminishes during late merger and vanishes altogether following merger. By contrast, the viscous 
torque drives the flow into the strong-field region and onto the primary during all epochs. Following 
binary coalescence, the viscous torque alone governs the time-dependent accretion onto the remnant, 
as well as the temporal behavior, strength and spectrum of the aftermath electromagnetic radiation 
from the disk. We solve our GR-hybrid equation for a representative BHBH-disk system, identify 
several observable E&M signatures of the merger, and compare results obtained for the gas and 
E&M radiation with those found with the Newtonian prescription. 

PACS numbers: 98.62.Mw, 98.62.Qz 



I. INTRODUCTION 

Binary black hole (BHBH) mergers are likely to occur 
in regions immersed in gas, and the capture and accretion 
of the gas by the binary may result in appreciable elec- 
tromagnetic radiation. There exists the realistic possi- 
bility of detecting electromagnetic "precursor" radiation 
prior to the merger and before the maximum gravita- 
tional wave (GW) emission from a BHBH merger [i|, y|- 
Then, follovifing the detection of gravitational waves, ob- 
serving electromagnetic "afterglow" radiation could pro- 
vide further confirmation of the coalescence [3|-[9| . Such 
electromagnetic radiation can also serve as a useful probe 
of the gas in galaxy cores or in other regions where merg- 
ers take place, as well as a diagnostic of the physics of 
black hole accretion. This diagnostic may be particularly 
revealing once the masses and spins of the merging com- 
panions and BH remnant are determined from the GW 
signal. 

In this paper we focus on a geometrically thin, Keple- 
rian disk orbiting in the plane of a spinning BH with 
a low-mass companion. To follow the orbit-averaged, 
secular evolution of such a BHBH-disk system, a sim- 
plified, vertically integrated, 1 -I- 1 - dimensional Newto- 
nian model equation [Eq. ([T]) belowlhas been adopted 
in many previous studies (see, e.g., [l|, 0, [l3 and refer- 
ences therein). We have demonstrated how the steady- 
state solution to this equation can be used to determine 
the disk structure and electromagnetic radiation spec- 
trum during the long inspiral epoch prior to binary-disk 
decoupling, during which a quasistationary treatment is 
applicable ^11]. To illustrate this approach we solved the 
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steady-state equation for representative BHBH-disk sys- 
tems at decoupling, employing simplified prescriptions 
for the required viscosity v{r) and disk scale-height h(r ) 
profiles. Our steady-state approach was extended in [l2| . 
where the Shakura-Sunyaev [13] "one-zone" prescription 
for radiation transport was adopted in conjunction with 
a "^"-disk viscosity law to obtain these profiles self- 
consistently (see also [ij]). 

Here we provide an alternative evolution equation that 
better approximates the strong-field, relativistic nature 
of the circumbinary disk. This equation treats the vis- 
cous torque in full GR for gas flow in a thin Keplerian 
disk. The disk, which is not self-gravitating for densities 
of interest here, evolves in the background spacetime de- 
termined by the more massive primary, assumed to be 
a (quasi-)stationary Kerr black hole. The tidal torque, 
arising from the presence of the low-mass secondary, is 
handled in the Newtonian limit. The later approximation 
is reasonable since the tidal torque is strongly peaked in 
and just outside a narrow gap in the disk centered on 
the orbit of the secondary. This torque plays its most 
important role prior to binary-disk decoupling, when the 
binary separation is large and lies outside the strong-field 
region of the primary. Moreover, the tidal torque disap- 
pears altogether following merger. By contrast, the vis- 
cous torque drives gas into the strong-field region during 
all epochs, including the post-merger phase. A GR treat- 
ment of the viscous-driven accretion onto the primary 
and the post-merger remnant is particularly important 
for making predictions of any observable, 'precursor' and 
'aftermath' electromagnetic radiation that may accom- 
pany the GW burst. 

A fully reliable description of the accretion flow and 
associated radiation really requires a radiation magneto- 
hydrodynamics (MHD) simulation in full general relativ- 
ity in the 3 + 1-dimensional, dynamical spacetime of the 



merging BHBH binary. Newtonian hydrodynamic simu- 
lations incorporating some of the relevant physics have 
been performed in various dimensions and levels of ap- 
proximation (see, e.g., P. la. Pzl. [l5l - [l7| . while GR simula- 
tions are in their preliminary stages (e.g., |18l426| ). Only 
recently have the first relativistic MHD simulations of 
a BHBH-disk system been performed: Noble et al. [27i| 
adopt post-Newtonian gravitation to perform simulations 
of an equal-mass system, excising the region inside the 
binary orbit, while Farris et al. |28| summarize simula- 
tions in full GR that cover the complete spatial domain, 
including the black holes. Both of these relativistic MHD 
simulations deal primarily with geometrically thick (i.e. 
warm) disks in which an "effective viscosity" is provided 
by magnetic fields driven turbulent by the magnetorota- 
tional instability (MRI). 

The model discussed here is mainly relevant for geo- 
metrically thin (i.e. cool) circumbinary disks. Although 
based on a simplified orbit-average description, it should 
delineate many of the qualitative features characteriz- 
ing the evolution of BHBH-thin disk systems. Also, the 
model may be useful for selecting input parameters and 
identifying scaling behavior for future, more detailed nu- 
merical simulations. In addition, the resulting solutions 
can provide approximate initial disk profiles for such sim- 
ulations. It is in this spirit and toward these purposes 
that we propose the adoption of our simple GR-hybrid 
equation. We hope that it provides a starting point for 
improved GR modeling along these lines. 

In Section II we review the Newtonian binary-disk 
model and the required elements that enter the secular 
evolution equation. We also summarize how the result- 
ing accretion rate onto the primary and the local elec- 
tromagnetic flux and total luminosity from the disk can 
be calculated. Simplications that arise when describing 
the pre-dcoupling and post-merger epochs are summa- 
rized. In Section III we present the GR-hybrid model 
and retrace our previous discussion, now adapting it to 
the GR-hybrid equation. In Section IV we provide a nu- 
merical example by solving the equations for an illustra- 
tive BHBH-disk system. We begin our integrations prior 
to binary-disk decoupling and proceed through inspiral 
and merger, comparing the Newtonian and GR-hybrid 
solutions. In Section V we outline future work that will 
improve the model. We adopt geometrized units and set 
G = 1 = c throughout. 
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Here G = — Tvis + Ttid is the total torque, Tvis is the vis- 
cous torque, Ttid is the tidal torque on the disk from the 
presence of the secondary, and Q = ^k = (M/r^)^/^ is 
the Keplerian orbital frequency about the primary, cen- 
tered at r = 0. The mass of the primary is M and the 
secondary qM, where q <^ 1. The viscous torque density 
is given by the standard equation [l|, y, [2^, [SO] 



dr dr \ dr 



(2) 



We approximate the (orbit-averaged) tidal torque den- 
sity by using the expression adopted by Armitage and 
Natarajan [1| 



dTtid 
dr 

where A(r, a) is given by 
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+ {fq'M/2r) (a/ Apt , r>a 
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In Eq. ([4]) / is a dimensionless normalization factor and 
Ap is given by Ap ~ maxdr— a|, h), and a(t) is the orbital 
radius of the secondary. Calibrating the above expression 
for the tidal field against high-resolution, hydrodynami- 
cal simulations in two-dimensions for a low-mass, black 
hole secondary interacting with an outer accretion disk, 
Armitage and Natarajan find that the value / ~ 0.01 
best fits the simulation results. Equations ([3]) and ^ fur- 
nish a reasonable analytic approximation to the results 
obtained from summing over the pointlike contributions 
from the Lindblad resonances in the disk [31, ^J . (Simi- 
lar, but slightly different, forms for the tidal torque also 
have been used in the literature; see, e.g., |2l. I33l435| . 

Assembling the above expressions then yields the final 
evolution equation 
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The accretion rate through any radius r in the disk 
may be calculated from 



II. THE NEWTONIAN EVOLUTION 
EQUATIONS 

A. Disk Evolution 

For reference and comparison we write down the stan- 
dard Newtonian evolution equation for the surface den- 



M{t,r) = 27rrE(-Wr) 

d{r^n) 



dr 

Combining Eqs. ([T]) and ^ yields 
5S 1 d 
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To solve the above evolution equation we impose the 
following boundary conditions; 



b.c.'s 
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In Eq. ([8|) rout is the outer radius of the disk and rjsco is 
the ISCO radius of the primary. Typically, rout ^ risco 
and in some cases one can take rout — >■ oo. We retain 
the solution for finite rout in part to facilitate numerical 
implementation of the outer boundary condition. 



B. Orbital Evolution 

The rate at which the secondary black hole migrates 
inward is determined both by back-reaction to the tidal 
torquing of the disk and by gravitational wave emis- 
sion [i,E3i: 



da/dt — {da/dt)tid + (da/dt)Q^N 



(9) 



where 



W^ih.--jjm-,lj-'^^'"- (-7;;;;). CI) 

and where (da/dt) qw is given by the familiar 
quadrupole-radiation orbital decay law, 



[da dt) Gw = — = ^ 

5 w^ 



tow 



(11) 



where C = 4q/(l -I- q)"^. In Eq. (|10l) the integration is 
over the entire disk, although most of the contribution 
from tidal torques arises close to the gap boundaries near 
r fa a. 

During any epoch in which back-reaction to tidal 
torques is not important, Eq. ^ can be integrated to 
yield 

a{t)/a{0) = (1 - 4t/tGw(0))'/^ , icwAtid « 1 , (12) 
where t = marks the beginning of such an epoch. 

C. Electromagnetic Radiation 

The local radiated emission from the disk arises both 
from viscous and tidal dissipation. We assume that all 
of the dissipation is radiated locally, whereby the local 
electromagnetic flux F{t, r) from each side of the disk is 
equal to the local dissipation rate D(t, r) per unit surface 
area. The rate of viscous dissipation is (29| 

9 M 
D^Ut, r) = -1.1]— = Fvi.(t, r). (13) 

The rate of tidal dissipation is given by [Ifl, |3g| 

Aid(i, r) = i ma) - nir)) AS = i^tid(i, r). (14) 



The local flux generates the luminosity L{t, r) accord- 
ing to 



Fvis(i,r) = M2i^vis(i,r)/Afeq 

1 /mV d 
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(Lvi,(i,r)/Meq), (15) 



with a similar relation between Ftid{t,r) and itid(^, *")• 
The total local flux is then given by 



F{t,r)^F^i,{t,r)+Ftid{t,r), 



(16) 



and the total luminosity integrated over the entire disk 
(both sides) is 



L(i)=ivis(0+itid(t). 



(17) 



In Eq. (jl5|) Afcq is the accretion rate in an equilibrium 
disk about a single black hole of mass M. When we 
compare Newtonian and GR-hybrid results we will use 
Eq. ([55)1 for Mcq in all of our normalizations. 



D. Quasistationary Solution: Pre-Decoupling 

As the binary inspiral proceeds from large separation, 
the inspiral timescale due to gravitational wave emission 
eventually becomes shorter than the viscous timescale 
in the disk, at which time the binary decouples from the 
disk and ultimately merges. We define the decoupling ra- 
dius ad to be the separation at which the two timescales 
become equal. Prior to BHBH-disk decoupling the bal- 
ance between tidal and viscous torques drives the disk 
to a quasistationary equilibrium state, perturbed slightly 
by small amplitude, spiral density waves emanating from 
the edges of the gap. Previously we solved the disk evo- 
lution equations in steady state to determine the quasis- 
tationary, (orbit-averaged) surface density profile prior 
to decoupling as a function of the the binary separa- 
tion [11[; see also [13, 'l4!|- For these early epochs we 
set a = constant and dYi/dt = in Eq. ([S]) to obtain 
the density profile. This quasistationary solution is used 
below as initial data for the Newtonian time-dependent 
simulations that evolve the binary-disk system from pre- 
to post-decoupling, continuing all the way through the 
late inspiral, merger and post-merger phases. 

The accretion rate M in steady state is independent of 
r. In steady state Eq. ([5]) admits a first integral which, 
when combined with Eq. ([6]), yields a first-order ODE, 



Af = 27r 



3r 



1/2 



d(r^/^vY) 2AJ:r^/'^ 



dr 



M^l^ 



constant. 



(18) 

We could solve the second-order elliptic equation ob- 
tained by setting the right-hand side of Eq. ([5]) to zero 
to obtain the steady-state density profile, then evaluate 
Eq. (IT51) for the accretion rate. Alternatively, we could 
integrate the first-order Eq. (|18p directly for the density. 



which, when the boundary conditions are implemented, 
automatically provides M as an eigenvalue. We chose 
the later strategy in [ll| , but adopt the former approach 
in Section IIVB II in obtaining the initial data. 



III. THE GR-HYBRID EVOLUTION 
EQUATIONS 

A. Disk Evolution 



E. Quasistationary Solution: Post-Merger 

Following binary merger the tidal torque vanishes while 
gas in the disk continues to diffuse inward on a viscous 
timescale, accreting onto the remnant black hole and ulti- 
mately settling into a final, stationary equilibrium state. 
This stationary disk configuration is described by well- 
known analytic density and temperature profiles, as well 
as analytic local fluxes and distant total luminosities, and 
these quantities provide useful checks on the late stages 
of any disk evolution calculation. The final equilibrium 
density profile, obtained by integrating Eq. ([5]) in steady- 
state in the absence of the tidal torque, is given by the fa- 
miliar result for a Shakura-Sunyaev Newtonian thin disk 
around a single black hole (see, e.g., [30| and references 
therein) , generalized for a disk of finite radial extent [11| : 
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The second equality above thus yields the steady-state 
accretion rate Meq in terms of the density and viscosity 
at the outer boundary: 



Afoq = STTi^outSo 



1 



' , _ 1/2 , 1/2 ■ 
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(20) 



The corresponding stationary fiux, due entirely to viscous 
dissipation, may be expressed as 



^vis(.) = ^f^) {y-^Jr^l^. (21) 



Stt V r 



The flux, together with Eq. ([TSl) . gives the differential 
luminosity. 



(ivis(r)/Afe, 



.,„,..-- V 4f (>--->"^)- <^^' 



The total steady-state luminosity integrated over the en- 
tire disk is then given by 



-'^vis — 



1-3^^ + 2- 
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(23) 



We propose the following GR-hybrid evolution equa- 
tion for the rest-mass surface density (E = J p^dz, where 
po is the rest-mass density) to replace Eq. ([5|): 



9E 
'dt 



Tr dr 



Q^ dr 



■i/2^s: 



2AS]r3/2 
Ml/2 



(24) 

In assembling the above equation we specialized to the 
Kerr metric in Boyer-Lindquist coordinates to describe 
the (quasi-) stationary spacetime established by the more 
massive primary black hole. We used this metric to ex- 
press the following functions, many of which were intro- 
duced by Novikov and Thorne [33] (see also Page and 
Thorne Q): 

M = mass of primary black hole, 
J = spin angular momentum of the primary hole. 



J/M^ 



< a* < 1, 



a; = (r/M)i/2, 

r = ^^-i/2. 



= 1 + a, a; 
= l-3a;"^ 



2a*a; 



-3 



S)^l-2x-^ +alx-^, 
^ = 1 ~ 2x^'^ + a^x^^ , 
^ = Eq. (35) in |38|, 



(25) 



We note that in the case of a thin disk around a sta- 
tionary Kerr black hole E as defined above is a scalar 
invariant, like po. 

Our proposed disk evolution Eq. (|24p has the following 
features: 

1. We assume that during all epochs the gas fiow takes 
place in the background geometry of the more mas- 
sive primary, which we approximate by the station- 
ary Kerr metric. The viscous torque, described by 
the first term on the right-hand side, is treated in 
full GR for gas flow in a thin Keplerian disk. In 
the absence of the tidal torque term arising from 
the presence of the secondary, the equation reduces 
identically to the evolution equation presented in 
Lightman and Eardley [33 for time-dependent disk 
accretion onto a single Kerr black hole. 

2. The tidal torque density, accounted for by the sec- 
ond term on the right-hand side, is based on the 
Newtonian formula, Eq. ( ^ . 

3. The entire equation reduces identically to Eq. ([5]) 
in the weak-fleld region, i.e., for r/M 3> 1, whereby 
T, &,'^ and Q all appproach unity. 



4. In the absence of tidal torques, the steady-state 
solution of Eq. ([^ gives the same profile for the 
product i/T, that characterizes a standard relativis- 
tic Novikov-Thorne accretion disk around a single 
Kerr black hole [again generalized for a disk of finite 
radial extent; see Eq. (|32|)]. 



Basing the tidal torque on the Newtonian expression is 
motivated by the fact that this torque is strongly peaked 
near the orbit of the secondary and plays its most impor- 
tant role when the binary separation a is large. Specifi- 
cally, the disk radii in which the torque is most significant 
typically satisfy r ^ a 3> Af and thus reside in the weak- 
field region outside the primary. (We neglect any accre- 
tion onto the secondary, which is expected to be small). 
By contrast, the viscous torque drives gas into the strong- 
field region and into the primary during all epochs and 
this flow requires a full GR treatment Moreover, follow- 
ing merger, the tidal term vanishes and Eq. ([M)) reliably 
accounts for the inward diffusion of gas, the filling of any 
pre-merger gaps in the disk, the time-varying accretion 
onto the remnant, and the relaxation of the disk and ac- 
cretion rate to a (quasi-)stationary state, all in full GR. 

A more rigorous treatment would incorporate a rela- 
tivistic tidal torque density dTtid/dr in place of the New- 
tonian expression used here. Such a relativistic formula 
presumably can be obtained by employing the relativistic 
tidal torque derived by Hirata [40, |4l[ for a single Lind- 
blad resonance in an accretion disk that orbits a Kerr 
black hole and is perturbed by a small secondary. This 
formula may be summed over many resonances, treated 
as a continuum, to get a smooth torque density. Such a 
sum has been carried out only for a Newtonian disk [31j , 
and has been used here. However, as described above, 
employing this Newtonian formula in a first approxima- 
tion should be adequate to treat many of the epochs of 
interest during the merger event. 

To solve Eq. (P^ we impose the same boundary condi- 
tions as specified by Eq. ([8]) . The radius rjsco is given by 
the familiar expressions for a Kerr black hole in Boyer- 
Lindquist coordinates (see, e.g. [421, Eq. 12.7.24). 

The rest-mass accretion rate through any radius r is 
given by 

Mo{t, r) = 27rrE(-wO^^^^ (26) 
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Q" dr 
Combining Eqs. ([Ml) and (^5)) then yields 
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(27) 



27rrr dr 
26]) and ^ reduce to Eqs. dH) and in the New- 
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tonian limit. 



relativity, even in vacuum. Post-Newtonian (PN) ap- 
proaches, based on expansions in u^/c^, can treat most 
of the inspiral epochs, but break down once the orbital 
separation shrinks to within a few times the radius of the 
massive primary. Treating the problem without approxi- 
mation using the tools of numerical relativity is not com- 
putationally practical for following the inspiral from large 
separations characterizing binary-disk decoupling, but it 
can match onto PN trajectories at late times to continue 
the late-inspiral motion through plunge, merger and ring- 
down (for an overview and references, see (43|). However, 
numerical relativity cannot yet evolve binaries with mass 
ratios q < 10~^ because of the excessive dynamic range 
and resulting resolution requirements. For such small 
mass ratios black hole perturbation theory provides the 
best approximation, although it is computationally and 
analytically expensive. One approach involves the calcu- 
lation of the self-field acting on a test particle and fol- 
lowing how it alters the orbital trajectory (For a status 
report and references see Q). A simpler, but more ap- 
proximate, method is the "radiative-adiabatic" scheme, 
where the inspiral is treated as a sequence of adiabati- 
cally shrinking geodesies. The shrinkage is determined 
by effectively calculating the rate of change of the con- 
stants of the motion (energy, angular momentum and 
Carter constant) due to GW emission. (For a summary 
and references see [4^ ) . Simplifications arise for the case 
of interest here, where the orbit is nearly circular and 
resides in the equatorial plane of the primary. 

It is therefore possible to modify Eq. ([TT|) to obtain 
a more reliable expression for the orbital decay due to 
GW emission that accounts for higher-order general rel- 
ativitistic effects. Here, however, we will continue to use 
Eq, ([TT|) for simplicity, as it is adequate to illustrate disk 
evolution via the GR-hybrid approach in a first approx- 
imation and can be generalized in subsequent analyses 
using the GR methods summarized above. Moreover, 
the deviations from the simple quadrupole approxima- 
tion that arise when the secondary approaches the pri- 
mary only lasts for a brief time interval, during which the 
bulk of the disk barely alters its structure. For the same 
reasons we will continue to use Eq. ()IOp to approximate 
the back-reaction of the tidal torque on the companion. 



C. Electromagnetic Radiation 



The local radiation flux F^g"" 



B. Orbital Evolution 



^{t,r) removes the local 
viscous dissipation in the disk. Measured from each side 
of the disk per unit surface area by an observer comoving 
with the gas, it is given by 



The inspiral of a low-mass black hole companion onto a 
more massive primary is a nontrivial problem in general 






1/2 



5r(i,r)|, (28) 



where Wit, r) is the vertically-integrated shear 
stress [37| . 

/■ 3 M^/^ 'Si 

W{t,r)^Jdzt^, = -.i:-^-. (29) 

We then define 



^K^)7lV.... 



/Mcq, (30) 



where Meq is the rest-mass accretion rate in an equi- 
librium disk about a single black hole of mass M (see 
Eq. I33p . Consistent with our adopting Newtonian ap- 
proximation for the tidal torque and tidal dissipation, 
we again use Eq. (|T4l) for i^tid and define Ftid{t,r) by 
analogy with Fvis{t,r). 

Following [4a] we define a flux F in terms of the 
radiated energy as measured by a distant, stationary 
observer: F{t,r) = dE / {rdrd(j)dt) = -utF''°"'{t,r). 
Here E is the energy measured by the distant observer 
and u" is the 4-velocity of a fluid clement in a circu- 
lar equatorial geodesic orbit about the primary: ut = 
—u)° ■ d/dt = — ^Z*^^/^. Now the spacetime metric, 
which is dominated by the primary, is stationary. More- 
over the orbit-averaged disk and its associated electro- 
magnetic emission evolve on a slow, secular timescale 
(~ min[ivisjiGw]) during most phases. In this limit the 
total luminosity measured by a distant observer may be 
computed from L{t) = 2dE/dt = Att f F{t' ,r'ydr' or 
d(L/Moq)/d(lnr) ^ Anr'^ F {f , r) / M^q, where t' is the re- 
tarded time from the observer to the source (the "fast- 
light" approximation). We omit the small correction for 
any emitted radiation captured by the black holes. 



D. Quasistationary Solution: Pre-Decoupling 



The discussion in Section IIIDI again applies. For bi- 
nary separations a ^ a^, the disk profile relaxes ap- 
proximately to a quasistationary profile found by setting 
dYj/dt = in Eq. p4)) and solving the resulting elliptic 
equation. We do so below in Section HVB II to determine 
the initial data for an evolution calculation. In this limit, 
the steady-state accretion rate Mq given by Eq. (pS)) sat- 
isfies 



A'/o = constant. 



E. Quasistationary Solution: Post-Merger 

Following black hole merger, a transient epoch ensues, 
wherein the gas diffuses inward toward the primary ac- 
cording to Eq. (|24)) in the absence of tidal torques and fills 
in any gaps that had been generated by the secondary. 
The disk eventually settles into a steady-state, relativis- 
tic Novikov-Thorne thin disk about the remnant black 



hole, whereby the surface density satisfies dT,/dt — 0, 
yielding 



r^3/2 ^2 



Mcq-^^/^ 



^out 



[post — merger] 



- 3. ^2 ^^ (32) 

Note that while Eq. ((32|) reduces to the Newtonian re- 
sult, Eq. (I19p . as r — >■ oo, the equilibrium profiles differ to 
C'(M/r)^/2^ a significant difference for gas near the rem- 
nant black hole. The second equality in Eq. ([5^ gives 
the steady-state, rest-mass accretion rate Meq in terms 
of the density and viscosity at the outer boundary: 



Meq = StTI^o 



^l 



"avA. "^out 



(33) 



The gas moves in a nearly circular geodesic orbit with 
an angular velocity 



^ 



r3/2 ^' 



and an inward radial drift 



3/2 



(34) 



(35) 



as measured in the orthonormal orbiting frame. Com- 
bining Eqs. (pO]) and ([32]) yields the comoving stationary 
fiux, which is now due entirely to viscous dissipation: 



^ VIS in - ^^\ ^ \ ,^1/2- 



(36) 



in \ r 

The fiux results in a steady-state, differential luminosity, 
d ( ^ ,.,,-. \ 3M^^ 



dlni 



(Lvis(r)/Me, 



2 r 



(37) 



Integrating Eq. ((37)) over the entire disk yields the total 
observed luminosity Lvis- For an infinite disk this inte- 
gration yields 



-^vis/-''^^cq — t -C/isco 



(= V)^ ''out -^ oo. 



(38) 



where i^isco is the binding energy per unit mass of a test 
particle in a circular geodesic orbit at risco, 



-^isro — 



rfsco - 2Mn, 



a^M^/Wn, 



>(4co - SMrisco + 2a*MyMf~;)i/2 



(39) 



(31) For a large, but finite disk with M <C Tout < oo we have 



Lvis/ Ml 



cq 



1-^i. 



l-Eu 



DvisZTrr dr 



3 M 
2 r-out 



(40) 



The right-hand side of Eq. ([55)1 yields the well-known 
efficiency 77 of stationary accretion from an infinite disk 
onto a Kerr black hole: 5.72% for a* = and 42.3% for 
a.^, — 1. The efficiency is less for a finite disk, as indicated 
by Eq. ^. 



IV. NUMERICAL EVOLUTION 



A. Nondimensionalization 



To solve Eq. (|24)) numerically it is convenient to intro- 
duce the same nondimensional variables defined in 11 ll : 



s = (r/rout)^^^, si = (a/rout)^^ J S2 = (nsco/»'out)^^^, 
E = E/Sout, v = v/i'out, y = sT,, 
h = h/r, T = t/2Uisirout)- 



Here 



tvis{r) = - — 



(41) 



(42) 



is the characteristic viscous timescale at radius r in the 
disk. 

In terms of these variables, equation (l24l) becomes 



dy Id 
dr Fs^ ds 

where 






vy- 



9*{s)y 



9*{^) = 



gsl 



s > Si 



-gs s < si 



(43) 



(44) 



and where 



2 , 2 /^out\"^ /'rout\ 



1/2 



(45) 



Eq. (|43|) must be solved for s E [s2, 1] subject to the 
boundary conditions 



b.c.'s 



y 



1, .s = 1 

0, S = S2 



(46) 



We integrate Eq. (j43p numerically, implementing a 
second-order, finite-difference, Crank-Nicholson scheme. 
Such an approach allows for arbitrarily large time-steps 
without the restriction of a Courant condition to insure 
stability. We choose a logarithmically increasing grid in 
radius to cover the large dynamic range in the disk with 
adequate spatial resolution everywhere. 



separation Ud by the condition iGw(ad) = /?ivis(2ad) (set- 
ting /? = 0.1), which yields 



Ud/M : 



128 
15-2' 



■^^ \ M J \ M J 



l/(n+2) 



(47) 



We next fix /i = h/r = 0.1, which essentially establishes 
the disk thickness near r — a, where it most matters. 
To set the scale for the density and disk size in physi- 
cal units we fix Sout and Tout, which determine the disk 
mass Mdisk- Finally, we set ^'out by specifying the final 
accretion rate onto the black hole remnant for an infinite 
disk, Mrem = STTfoutSout (sec Eqs. [201 or [32] in the hmit 
Tout 3> M) to be a fraction 7 of the Eddington value, 7 = 

M.cm/MEdd- Here MEdd = L^dd/v = 47rMTOp/(77crT), 
where rup is the proton mass, ax is the Thomson cross- 
section and 77 is the radiative efficiency [see Eq. 
This condition yields 



t/M 



47 mp 
3 77 or So 



(48) 



As in [il|, we assign the values n = 0.5 and q — 
5 X 10-3 and set / = 0.01, rout = lO^M, 7 = 0.1, and 
Scut = 1-5 X lO^gcm"^. Our choice of parameters gives 
Mdisk/Mo - 5 X 10^ M| (Ms = M/10^Mq), which is 
safely smaller than the mass of the secondary BH for 
all M < 2 X IO^^Mq. These values for the asymptotic 
disk density and disk mass are comparable to cases con- 
sidered in [1, 2]. We also set J/NP = 0.5, which gives 
Hsco = 4.23M. The value of the primary mass M scales 
out of the problem when solved in dimensionless form ac- 
cording to Eq. (|43|. 

The adopted parameters give a decoupling radius 
ad/M = 18.1 and a dimensionless tidal torque param- 
eter g = 0.0194 [see [U], Eq. (35), for definition]. 

A more sophisticated treatment could adopt a one-zone 
approach as in a Shakura-Sunyaev or Novikov-Thorne 
disk. There one employs a local radiation prescription 
that yields the local temperature and pressure for an 
a-disk or /3-disk viscosity law, and uses these to de- 
rive the viscosity and h/r profiles self-consistently (see, 
e.g. [IJ, [iJI). However, the simplified, but physically 
plausible, assignments chosen here are sufficient for illus- 
trating the implementation of the GR-hybrid equation 
(where we merely take it out for a "test-drive") and we 
postpone a more detailed analysis for a future investiga- 
tion. 



B. A Numerical Example 

To illustrate how a BHBH-disk system evolves when 
governed by the GR-hybrid equation we track a typi- 
cal BHBH-disk system by integrating this equation in 
time. Following our approach in [1 1| . we take the vis- 
cosity to have a power-law profile iy{r) ex r". We 
then specify the system by first choosing the parameters 
q, a/M, roat/M ^ 1 and n. We determine the decoupling 



1. Initial Data 

We start the evolution when the binary separation is 
at a/M — bod/M — 90.3. We thus begin to track the 
inspiral before binary-disk decoupling, when for much of 
the disk the evolution is still quasistationary. We can 
therefore set the initial density profile E to the quasista- 
tionary profile found by setting dYj/dt = in Eq. 



(i.e., dy/dr = in Eq. |43| and solving the resulting el- 
liptic equation, as discussed in Section UlID I 

At t = we have g = 9.68, where the quantity g is 
defined by Eq. (55) of \IV\ and is related to g according 
to 



1 



9 



>5 



a 

''out 



(l/2-«) 



h{a) 



(49) 



As shown in that reference, g measures the ratio of the 
tidal to viscous torque on the disk at r ~ a. Whenever 
g ^ a few prior to merger the dimensionsionless accretion 
rate satisfies m ^ 1, i.e., the accretion rate onto to the 
primary is much reduced below the value it would have in 
the absence of tidal torques from the secondary. At i = 
we find m = 1.03x 10~^, consistent with this expectation. 
Note that g decreases rapidly with increasing h. Hence 
for thicker disks with h ^ 0.1, g would be smaller and 
there could be little or no suppression of accretion by the 
secondary (cf. [12]). 

At i = we find icw/itid = 0.0588M8. This ra- 
tio, which measures the relative importance of tidal to 
gravitational radiation back-reaction forces on the sec- 
ondary, decreases as the inspiral proceeds. Hence tidal 
back-reaction is not important in this particular scenario 
for all M < 10^. Accordingly, we use Eq. (IT^ to track 
the orbital separation, consistent with the discussion in 
Section HHBl 

To assist in the interpretation of the numerical results 
and figures below, we list the following conversion of 
nondimensional to physical units applicable to this sce- 
nario: 
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FIG. 1. Snapshots of the disk surface density profile at 
selected times. Profiles are shown for the initial disk at 
r = {dotted red line), the final equilibrium disk at r = oo 
{solid black line) and for several intermediate times, r = 
0.03707, 0.04943, 0.05251, 0.05366 (pre-merger dashed black 
lines) and 0.05406, 0.06178, 0.09267, 0.2471 (post-merger dot- 
dashed blue lines). Decoupling occurs at at r = 0.05360 
{t = 4.099 X lO^Mg yrs), and merger occurs at r = 0.05368 
{t = 4.106 X 10^ Ms yrs). 



i(yrs) = 0.7649 x lO^Mgr, 
ad{an) = 17.84M8, 



Meq(M0 yr-i) 
Leq(erg s"^) 



= 7^Edd^out 
= 7-^Edd=^out = 



= 0.3028M8, 
1.409 X 10''^M8.(50) 



The subscript 'eq' refers to final quasistationary values 
associated with accretion onto the black hole remnant. 
The factor ^out = ■^out/('^^out -^out) corrects for a rela- 
tivistic disk which is finite and not infinite, as is the case 
here {cf. Eq. 



2. Evolution 

The surface density profile is plotted at selected times 
during the evolution in Fig. [TJ A gap forms in the disk 
near the orbital radius of the secondary at r = a{t) and 
moves with the secondary as it spirals inward. Disk- 
binary decoupling occurs after a{t) has reached a^, which 
happens at r = 0.05360 {t = 4.099 x lO^Afg yrs). Prior to 
that time, the density profile evolves in a quasistationary 
manner and remains close to the quasistationary solution 
obtained by setting dY,/dt = at each orbital separation 
a{t) > ad. Tidal torques, which are strongest near the 
orbital radius of the secondary, cause a pile-up of the 



matter and create a density peak just outside r — a{i). 
As the orbital radius shrinks, the peak surface density 
moves inward while steadily increasing. After merger, 
which occurs at r = 0.05368 (t = 4.106 x lO^Mg yrs or 
Ai = 6.57M8 yrs after decoupling), the tidal torques van- 
ish altogether while the residual viscous torques drive the 
inward diffusion of gas toward the remnant black hole. 
This inward diffusion reduces the peak value of the den- 
sity and eliminates the gap inside r = a^ altogether. By 
T = 0.2471 {At = 1.479 x lO^Mg yrs after merger), the 
density profile is seen in the figure to be nearly indistin- 
guishable from the equilibrium Novikov-Thorne solution 
for a relativistic disk accreting onto a single black hole, 
Eq. ((32| . The density at the ISCO vanishes at all times, 
as it is one of the boundary conditions. 

The rest-mass accretion rate as a function of radius is 
plotted at select times in Fig. [2j At i = t = the accre- 
tion rate, given by Eq. (1261) . is everywhere constant. Such 
a result is a consequence of demanding that d'S/dt — 
for the initial disk (see Eq. \T7]\ . The normalized value 
of the initial accretion rate, M/AI^q — 1.031 x 10~^, is 
much less than unity, the value characterizing an equil- 
brium disk with the same asymptotic density and viscos- 
ity about an isolated black hole. As discussed above and 
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FIG. 2. Snapshots of the rest-mass accretion rate at selected 
times. Profiles are plotted at the same times shown in Fig.[T] 



ll| . such a suppression of accretion flow is expected 



due to the high initial value of the torque parameter g. 
As the inspiral proceeds, the exact cancellation of the vis- 
cous and tidal torque terms in Eq. ([25]) breaks down and 
the mass flux grows behind the secondary. Soon after 
decoupling, the flow rate reaches values that actually ex- 
ceed the final equilibrium value for an isolated black hole 
by almost a factor ~ 10 near r ^ a^, as the tidal torques 
holding back the flow diminish in strength and the den- 
sity pile- up abates (the "bursting of the dam" ) . Follow- 
ing merger the accretion rate eventually settles down to 
the constant equilibrium value, Eq. ([55)1 . throughout the 
disk. 

It is interesting to compare the surface density pro- 
files determined by integrating the GR-hybrid and 
Newtonian evolution equations for the same primary 
and secondary black hole masses, disk parameters (i.e. 
Sout,;^out,nsco,?'out,'^, / and h/r) and secondary orbit 
a{t). In each case the initial data is determined by set- 
ting dY,/dt = in their respective evolution equation 
and solving the resulting elliptic equation. The compar- 
ison is provided in Fig. [3] The initial quasistationary 
profiles are nearly identical, except in the strong-field re- 
gion r/M < 20, where the Newtonian profile is higher. 
(A comparison of profiles in the strong-field region is 
of course influenced by gauge effects arising from the 
choice of radial coordinate, but E is a scalar invariant). 
The initial accretion rates are also slightly different (e.g. 
M/M^^ = 0.8892 X 10~^ in the Newtonian case), since 
the elliptic equations that determine the rates are differ- 
ent in the strong-field region. Prior to merger but after 



FIG. 3. Comparison of GR-hybrid and Newtonian disk sur- 
face density profiles at selected times. Profiles are shown for 
the initial disk at r = [dotted lines), the final equifibrium 
disk at r = oo [solid lines) and at one intermediate time r 
= 0.05366 (pre-merger dashed lines). GR-hybrid lines are in 
black and Newtonian fines in red; for each line type the lower 
(upper) curves are the GR-hybrid (Newtonian) lines. Merger 
occurs at r = 0.05368 [t = 4.106 x lO^Afg yrs). 



decoupling, the density profiles outside the orbital radius 
are close, but inside that radius they continue to depart. 
After merger, the equilibrium profiles in the strong field 
regime remain different: the peak value of the surface 
density is higher in the Newtonian case by 35%. The dif- 
ferences between the Newtonian and GR-hybrid solutions 
become more pronounced as the primary spin increases 
and a* -^ 1. 

Profiles of the nondimensional comoving flux emerging 
from each side of the disk are plotted at selected times 
in Fig. m The evolution of the flux is correlated with the 
evolution of the disk surface density. The peak flux oc- 
curs just outside the orbital radius of the secondary prior 
to merger, and the regions where the flux dips correspond 
to the density gaps near that orbit. Prior to merger, the 
tidal torque-driven density pile-up causes the peak flux 
to increase with increasing time and decreasing orbital 
radius. Following merger the tidal torque vanishes and 
the flux begins to decrease everywhere. Eventually the 
flux, now generated by viscous dissipation alone, settles 
into the equilibrium state corresponding to steady accre- 
tion of gas in a thin, relativistic disk onto a single black 
hole [Eq. §^] The flux vanishes at the ISGO. 

A comparison of the comoving fluxes determined from 
the GR-hybrid and Newtonian evolution equations is pre- 
sented in Fig. [51 Prior to merger the flux profiles are 
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FIG. 4. Snapshots of the comoving flux profile at selected 
times. Profiles are plotted at the same times shown in Fig.[T] 



FIG. 5. Comparison of GR-hybrid and Newtonian disk co- 
moving flux profiles at selected times. Profiles are plotted at 
the same times shown in Fig. [S] 



quite comparable but after merger the final equilibrium 
profiles to which the disks relax differ significantly in the 
strong-field region r/M ^ 15. The peak value of the final 
equilibrium comoving flux is a factor of 2.1 times larger 
in the Newtonian case. 

The ratio of the contribution of tidal heating to viscous 
heating to the comoving flux is plotted as a function of 
radius at selected times prior to merger in Fig. [6] Tidal 
heating dominates over viscous heating near r — a{t), 
but decreases rapidly both at smaller and larger radii. 
[Note that exactly at r = a{t) the tidal dissipation van- 
ishes in accord with Eq. (J14p . hence the sudden dip in 
the curves]. This sudden fall-off is anticipated because 
the tidal torque, due to the presence of the secondary, 
decreases rapdily with distance from the secondary. 

The contribution from various radii in the disk to the 
luminosity measured by a distant observer is plotted in 
Fig. [7] at selected times. The evolution of this quan- 
tity follows the general trends already found for the co- 
moving flux shown in Fig. |4l The same differential lu- 
minosity function (up to our normalization) has been 
plotted in [46] for relativistic, Novikov-Thorne thin disks 
undergoing steady-state accretion onto single Kerr black 
holes, where they are compared with three-dimensional 
GRMHD simulations of thin disks (h/r < 0.1) that relax 
to steady-state; see their Fig. 1. Our post-merger lumi- 
nosities are all driven to these Novikov-Thorne solutions 
at late times. The GRMHD results are in general agree- 
ment with these solutions, but do exhibit some emission 
inside the ISCO, plus a small inward shift of the peak 
emission to lower radii (see also [4711). These differences 
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FIG. 6. Snapshots of tidal-to-viscous comoving fiux profiles 
at selected times. Profiles are plotted at the same times shown 
in Fig. [T] prior to merger; following merger tidal dissipation 
vanishes. 
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FIG. 7. Snapshots of the differential luminosity profile, 
which shows the contribution from various radii in the disk to 
the distant luminosity, is plotted at selected times. Profiles 
are plotted at the same times shown in Fig. [T] 



are not deemed significant enough to affect the accuracy 
of, e.g., the continuum- fitting method, which employs the 
Novikov-Thorne model to estimate black hole spins. 

The evolution of the total electromagnetic luminosity 
from the disk measured by a distant observer is plotted in 
Fig. [51 There it is seen that the tidal dissipation is always 
less important overall than viscous dissipation and van- 
ishes altogether following merger. The luminosity rises 
sharply after merger, reaching a peak at t = 0.0542 and 
decaying slowly thereafter to its final, equilibrium value. 
The peak value is a full 3.08 times larger than its final 
equilibrium value given by Eq.|40l Leq = 0.0806Moq. The 
rapid decline of the tidal torques following decoupling 
and the sudden inward drift of matter from r ~ a^ is re- 
sponsible for the overshoot in luminosity. The full width 
at half-maximum of the luminosity curve is At — 0.03126 
{At = 2.39 X lO^Ms yrs). The rise, fall and asymptotic 
flattening of the luminosity curve may provide an elec- 
tromagnetic signature that a BHBH merger has occured 
in a circumbinary disk. 

Also shown in Fig. [5] are the Newtonian evolution 
curves for the same quantities. The results are quali- 
tatively similar, but the equilibrium accretion rates and 
radiation efficiencies are different from the GR-hybrid so- 
lution (e.g., Leq = 0.117Afcq, Or 45% higher, in the New- 
tonian case) and the peak luminosity overshoot is some- 
what smaller (about 15% lower in the Newtonian case). 

Fig. [S] shows that the evolution of the total electro- 
magnetic luminosity is correlated with the evolution of 



FIG. 8. Variation of the distant total disk luminosity with 
time. The dashed lines show the contribution from tidal dis- 
sipation, the dotted lines from viscous dissipation. The solid 
lines show the total luminosity. GR-hybrid lines are in black 
and Newtonian lines in red; for each line type the upper 
(lower) curves are the GR-hybrid (Newtonian) lines. Merger 
occurs at r = 0.05368 (i = 4.106 x 10^ Mg yrs). 



the accretion rate at the ISCO of the primary black hole. 
The overshoot of the accretion rate above the final equi- 
librium value, and its subsequent decay to the Novikov- 
Thorne value, drives the same time variation seen for the 
total luminosity. 

As in the case of thin-disk accretion onto a single, sta- 
tionary black hole, Newtonian and GR models for low- 
mass BHBH-disk systems typically give the same qualita- 
tive results for the observable E&M radiation. Newtonian 
calculations are thus sufficient to identify characteristic 
E&M luminosities and wavelengths. But the numerous 
factors of a ^ few that comprise the quantitative dif- 
ferences between the models are important if one hopes 
to use detailed observations to infer BH spins and other 
system parameters. For more complicated scenarios than 
the one analyzed here there can even be qualitative dif- 
ferences involving relativistic effects that must be taken 
into account (e.g., binary remnant recoil, misaligned BH 
spins, accretion jets, etc.). 



V. FUTURE WORK 

The numerical scenario summarized here is presented 
as a simple demonstration of the use of the GR-hybrid 
approach to track the orbit-averaged evolution of a thin. 
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FIG. 9. Variation of the distant total disk luminosity {solid 
lines) and ISCO accretion rate {dotted lines) with time. GR- 
hybrid lines are in black and Newtonian lines in red; for each 
line type the upper (lower) curves are the GR-hybrid (Newto- 
nian) lines. Merger occurs at r = 0.05368 {t = 4.106 x W^ Mg 
yrs). 



Keplerian disk orbiting a low-mass BHBH, accounting for 



some of the most important effects of general relativity. 
More detailed microphysics, combined with a parameter 
survey, will be necessary to fully explore the consequences 
of this model. Future applications should incorporate the 
following: 

1. A self-consistent treatment of the viscosity and h/r 
profiles by implementing a Shakura-Sunyaev-Novikov- 
Thorne one-zone description of each ring in the disk, 
employing a local radiation prescription together with 
an a-disk or /3-disk law for the viscosity to obtain the 
required profiles; 

2. A calculation of the observed radiation spectrum, 
by employing a ray-tracing or Monte-Carlo technique, 
adapted to a time-dependent relavistic disk in curved 
spacetime (see, e.g. |48l - [50| ). 

Several GR refinements can be implemented to im- 
prove the model while retaining the spirit of an orbit- 
averaged description of an evolving BHBH-thin disk sys- 
tem. They include the following: 

1. The replacement of the Newtonian formula with a 
fully relativistic expression for the tidal torque density, 
along the lines discussed in Section UlIAl 

2. Employing the true relativistic inspiral trajectory 
for the low-mass BH companion, calculated using one of 
the GR techniques discussed in Section UlIB I 

We intend to implement some of these improvements 
and apply the resulting formalism in future studies. 
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